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In  a previously-developed  phase-field  model  of  solidification  that  includes  convec- 
tion in  the  melt  [1],  the  two  phases  are  represented  as  viscous  liquids,  where  the 
putative  solid  phase  has  a viscosity  much  larger  than  the  liquid  phase.  In  this 
paper  we  report  numerical  computations  on  a simplified  form  of  this  model  which 
represents  the  growth  of  a two-dimensional  dendrite  in  a thin  gap  between  two 
parallel  thermally  insulting  plates.  In  these  computations  flow  in  the  liquid  arises 
because  of  the  differing  densities  of  the  solid  and  liquid  phases. 


1 Introduction 

The  notion  of  representing  the  boundary  between  two  bulk  thermodynamic  phases 
as  a diffuse  interface  dates  back  to  the  work  of  Poisson  [2],  Gibbs  [3],  Maxwell  [4], 
Rayleigh  [5],  van  der  Waals  [6]  and  Korteweg  [7]  in  the  19th  Century.  The  central 
assumption  is  that  there  is  an  interfacial  region  of  small  but  nonzero  thickness 
separating  the  two  bulk  phases.  In  such  models,  quantities  such  as  surface  tension, 
that  in  a sharp-interface  description  are  regarded  as  localized  on  the  interfacial 
surface,  are  instead  recognized  as  being  distributed  through  the  interfacial  region. 

Diffuse  interface  models  may  be  based  on  an  extended  thermodynamics  that 
incorporates  effects  involving  gradients  of  the  thermodynamic  variables  ( “nonclas- 
sical  terms”)  to  account  for  nonlocal  effects.  That  a model  incorporating  a diffuse 
interface  in  this  way  is  referred  to  as  “nonclassical”  is  perhaps  ironic  in  light  of  the 
above  history,  and  speaks  volumes  for  the  success  of  the  “classical”  sharp-interface 
description  of  interfacial  free  boundary  problems. 

Despite  the  overall  success  of  the  classical  approach,  there  are  still  special  sit- 
uations in  which  a diffuse-interface  description  of  an  interface  between  two  bulk 
phases  is  a viable  and  even  necessary  approach.  At  least  three  types  of  situations 
may  be  identified,  (i):  The  thickness  of  the  interface  becomes  comparable  to  or 
larger  than  other  mesoscopic  length  scales  of  interest  in  the  problem.  An  example  of 
such  a situation  is  in  the  case  of  a fluid  near  its  critical  point,  where  the  thickness  of 
the  interface  diverges.  Early  diffuse  interface  models  were  developed  to  investigate 
this  problem  (see,  e.g.,  van  der  Waals  [6]  and  Rowlinson  and  Widom  [8]).  (ii):  The 
length  scales  of  interest  in  the  problem  under  consideration  are  so  small  that  they 
are  comparable  to  the  thickness  of  the  interface.  Contact  line  problems  in  fluid 
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mechanics  (e.g.  Davis  [9])  are  potentially  in  this  category,  as  the  diffuse  nature  of 
a fluid-fluid  interface  may  become  important  at  the  small  scales  of  interest  near  a 
contact  line.  In  fact,  recent  calculations  using  diffuse-interface  descriptions  suggest 
that  the  force  singularity  associated  with  the  classical  free  boundary  description  of 
a moving  contact  line  (Dussan  V.  and  Davis  [10])  can  be  relieved  when  a nonzero 
interface  thickness  is  taken  into  account  (Jacqmin  [11,12],  Seppecher  [13]).  (iii): 
A diffuse-interface  formulation  becomes  a viable  computational  alternative  to  the 
classical  free  boundary  problem  when  the  morphology  of  the  interface  becomes 
very  complicated  or  changes  its  topology.  An  important  example  is  the  solidifica- 
tion of  dendrites,  where  sidearms  branch  from  the  main  stem  of  the  dendrite  in 
a complicated  dynamical  process  that  involves  both  growth  and  subsequent  coars- 
ening behavior.  Many  successful  computations  of  dendritic  growth  have  now  been 
performed  [14-21]. 

Diffuse-interface  theories  have  been  developed  and  applied  successfully  in  a wide 
range  of  other  physical  situations  as  well,  such  as  superconductivity  [22],  liquid 
crystals  [23],  spinodal  decomposition  [24,25],  ordering  transitions  in  alloys  [26-28], 
and  a variety  of  hydrodynamic  phenomena  [29] . 

Our  interest  here  concerns  a phase-field  model  that  accounts  for  both  solidifica- 
tion and  fluid  motion.  This  work  extends  the  phase-field  model  of  the  solidification 
of  a pure  material  that  was  first  proposed  by  Langer  [30,31]  and  subsequently  de- 
veloped by  a number  of  researchers  [32-37].  Phase-field  models  provide  an  example 
of  a diffuse- interface  model  in  which  an  order  parameter,  0,  is  postulated  whose 
value  indicates  the  phase  of  the  system  at  a particular  point  in  space  and  time 
(in  this  paper  0=1  and  0 = 0 denote  the  solid  and  liquid  phases,  respectively). 
Langer  represented  the  free  energy  of  a single-component  system  by  a gradient 
energy  functional  of  the  form 

^=^{je2|V0|2  + /(«,T)W,  (1) 

where  e is  the  gradient  energy  coefficient  and  T is  the  temperature.  The  free  energy 
density,  /(0,  T),  has  a double-well  structure  with  respect  to  0 in  which  the  two  local 
minima  correspond  to  the  solid  and  liquid  phases.  Langer  proposed  the  following 
governing  equations  for  the  phase  field  and  temperature: 


5JF_ 

6<f> 


(2) 


kv2T+L% 


(3) 


where  1/M  is  a positive  constant  termed  the  mobility,  c is  the  heat  capacity,  k is 
the  thermal  conductivity  and  L is  the  latent  heat  per  unit  volume  of  the  material. 
This  phase-field  formulation  replaces  the  free-boundary  problem  associated  with 
the  sharp-interface  model  of  an  interface  by  a coupled  pair  of  nonlinear  reaction 
diffusion  equations.  The  location  of  the  interface  is  represented  by  the  level  set 
0=1/2. 

An  early  attempt  to  include  fluid  motion  within  a phase-field  model  of  solidifica- 
tion is  due  to  Caginalp  and  Jones  [38,39].  They  appended  the  inviscid  momentum 
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equation  and  the  continuity  equation  to  the  phase- field  model,  but  did  not  address 
the  issues  of  momentum  balance  in  the  solid  and  capillary  contributions  to  the  stress 
tensor.  Diepers  et  al.  [40, 41]  have  employed  the  methodology  of  two-phase  fluid 
flow,  where  <j>  is  interpreted  as  a solid  fraction.  Their  model  is  used  to  study  coars- 
ening in  a binary  solid/liquid  mixture  with  and  without  fluid  flow.  Tonhardt  and 
Amberg  have  also  performed  two-dimensional  numerical  studies  using  adaptive  fi- 
nite elements.  They  study  the  effect  of  a shear  flow  on  dendritic  growth  morphology 
and  show  preferential  side-branching  on  the  upstream  side  of  the  dendrite  [42,43]. 

In  this  paper  we  briefly  describe  a recently-developed  phase-field  model  [1]  which 
allows  for  convection  in  the  liquid  phase.  This  model  has  three  notable  aspects: 
first,  we  represent  both  the  solid  and  liquid  phases  as  Newtonian  fluids  in  which 
the  viscosity  of  the  putative  solid  phase  is  specified  to  be  much  larger  than  that 
of  the  liquid  phase.  Second,  the  interface  is  ascribed  an  anisotropic  surface  energy. 
Third,  the  phase  transition  is  considered  to  be  first  order.  These  features  are  non- 
standard for  a model  which  treats  the  two  phases  as  Newtonian  fluids,  but  is  in 
keeping  with  our  intention  to  model  a solid-liquid  system.  We  note  that  in  many 
solidification  applications,  a fluid  model  is  used  for  the  thermodynamic  description 
of  the  solid  phase,  in  that  the  elastic  properties  of  the  solid  are  ignored.  In  order  to 
obtain  the  desired  viscosity  variation  between  the  phases,  the  viscosity  is  assumed 
to  depend  on  the  phase  field,  4>.  The  anisotropic  surface  energy  is  achieved  by 
employing  the  generalized  £- vector  formalism  [44].  Unlike  previous  diffuse  interface 
models,  which  incorporate  fluid  motion  coupled  to  a conserved  order  parameter 
description  [29],  we  adopt  a nonconserved  order  parameter,  4>,  in  line  with  our  aim 
of  directly  extending  conventional  phase-field  models  of  solidification  to  account  for 
convection.  This  has  the  advantage  that  we  may  treat  quasi-incompressible  systems 
[45],  in  which  the  density  field  is  taken  to  be  a prescribed  function  of  <f>. 

We  sketch  how  the  model  may  derived  in  the  setting  of  irreversible  thermody- 
namics. The  quasi-incompressibility  assumption  restricts  the  form  of  the  thermody- 
namic potentials  that  may  be  employed  [45].  The  model  comprises  the  compressible 
Navier-Stokes  equations  with  a modified  stress  tensor  that  includes  additional  terms 
related  to  gradients  of  <f>,  an  energy  equation,  and  a phase-field  equation  involving 
a material  time  derivative  of  <fi.  We  go  on  to  describe  computations  based  on  a 
simplified  form  of  this  phase-field  model.  In  particular,  we  study  a configuration 
in  which  a dendrite  grows  into  an  undercooled  liquid  between  two  thermally  insu- 
lating plates.  This  allows  us  to  avoid  directly  solving  the  generalized  compressible 
Navier-Stokes  equations  by  adopting  a Hele-Shaw  approximation.  The  densities  of 
the  solid  and  liquid  phases  are  allowed  to  differ,  and  we  study  numerically  the  effect 
of  the  density-induced  flow  on  the  growth  of  the  dendrite. 


2 The  Model 

We  consider  a non-isothermal  system  consisting  of  a pure  material  that  may  exist  in 
two  distinct  phases.  We  follow  the  standard  phase-field  methodology  and  introduce 
a phase-field  variable,  <j>(x,t),  whose  value  indicates  the  thermodynamic  phase  of 
the  system  as  a function  of  position,  x , and  time,  t.  A solid-liquid  interface  is 
represented  by  a thin  layer  in  which  the  phase  field  varies  rapidly  between  zero 


Interfaces  for  the  Twenty-First  Century 


3 


(liquid)  and  unity  (solid).  The  governing  equations  are  derived  by  following  the 
formalism  of  irreversible  thermodynamics  [37,46^8]  as  described  below. 


2.1  Governing  Equations 


We  assume  that  the  total  entropy,  S,  in  a material  volume,  Q(f),  of  the  system  is 
given  by 


<5 


-s 

JQ{t) 


ps  - ^e%T 2(V0) 


dV, 


(4) 


where  p is  the  density  and  s is  the  entropy  per  unit  mass.  The  first  term  in 
the  integrand,  ps,  is  the  classical  entropy  density  per  unit  volume  and  the  second 
is  a nonclassical  term  associated  with  spatial  gradients  of  the  phase  field.  Here 
the  gradient  entropy  coefficient  es  is  assumed  to  be  a constant  for  simplicity,  and 
T is  a homogeneous  function  of  degree  unity.  The  function  T allows  for  a general 
anisotropic  surface  energy  of  the  solid-liquid  interface  and  allows  the  Cahn-Hoffman 
vector  formalism  for  sharp  interfaces  [49,50]  to  be  generalized  and  extended  to 
diffuse  interface  models  [44,51].  An  isotropic  surface  energy  results  from  the  choice 
T(V0)  = |V0|. 

The  total  mass,  A4,  linear  momentum,  V,  and  internal  energy,  £,  associated 
with  the  material  volume  are  assumed  to  have  the  form 


(5) 

(6) 
(7) 


respectively.  Here  u is  the  velocity,  e is  the  internal  energy  density  per  unit  mass 
and  is  the  gradient  energy  coefficient,  which  is  assumed  to  be  constant.  The 
thermodynamic  relations 

n de 

de  = T ds  + dp  + —d(j),  (8) 

p-2  Q(p 


e = Ts-p/p  + p, 


(9) 


are  assumed  to  apply  locally,  where  p is  the  thermodynamic  pressure  and  p is  the 
chemical  potential  (or  Gibbs  free  energy  per  unit  mass). 

The  physical  balance  laws  for  mass,  linear  momentum,  and  internal  energy  are 
given  by 


— + [ 
dt  JsQ(t) 


H 

o 

(10) 

dV  f 

— — n • m dA, 

dt  JdQ(t) 

(11) 

qE-ndA—  / n-m-udA, 

JsQ(t) 

(12) 
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respectively,  where  n is  the  outward  unit  normal  to  5Q.(t),  m is  the  stress  tensor,  and 
qE  is  the  internal  energy  flux.  The  momentum  balance  (11)  requires  that  the  rate 
of  change  of  the  total  momentum  of  the  material  volume  results  from  forces  acting 
on  its  boundary  6Q(t)  (for  simplicity  we  neglect  body  forces  such  as  gravity;  their 
inclusion  is  straightforward).  The  energy  balance  (12)  equates  the  rate  of  change 
of  the  total  internal  energy  of  Cl(t)  plus  the  energy  flux  through  its  boundary  to 
the  rate  of  work  of  the  forces  at  its  boundary. 

In  addition,  the  entropy  balance  takes  the  form 

dS  f r 

— + / qs  • n dA  = / spToddV,  (13) 

at  JdQ(t)  Jn(t) 

where  qs  is  the  entropy  flux  and  sprod  is  the  local  rate  of  entropy  production.  The 
second  law  of  thermodynamics  is  then  expressed  by  the  requirement  that  spTOd  is 
non-negative. 

To  proceed  we  recast  the  conservation  laws  (10)-(13)  as  differential  equations. 
These  are  used  to  express  the  local  entropy  production  in  terms  of  the  fluxes  m,  qE, 
and  qs,  as  well  as  D(f>/Dt.  We  then  identify  forms  for  these  quantities  which  ensure 
that  the  local  entropy  production  is  non-negative.  The  fluxes  that  result  from  this 
procedure  involve  both  classical  contributions  and  non-classical  contributions  that 
depend  on  V0.  In  addition,  we  obtain  an  evolution  equation  for  the  phase  field.  The 
details  of  this  procedure  axe  given  in  Ref.  [1]  and  result  in  the  following  governing 
equations: 


Du 

p-—~  = V • m, 
y Dt 

M^  = 4(T)v-[rW)«l-p|, 

= v • [*VT]  + e%V  • [r(V*)|]  ^ + ms  : W, 


(14) 

(15) 

(16) 
(17) 


where  1/M  is  a mobility,  m is  the  stress  tensor  [see  equation  (23)],  ms  is  a modified 
stress  tensor  [see  equation  (24)],  eF  is  the  Helmholtz  gradient  energy  coefficient 
given  by  <rF(T)  = e\  + Te |,  g(T,p,<f>)  = e - Ts  + p/ p is  the  Gibbs  free  energy 
per  unit  mass,  and  £ is  the  generalized  ^-vector  [44]  whose  components  are  defined 
by  — dT{p)/dpj,  where  we  have  written  p = V0.  The  density  of  the  two  bulk 
phases  may  be  different,  and  we  will  assume  that  p depends  on  <j>  alone, 


P(0)  = Psr((t>)  + Pl[  1 - r(^)]. 


(18) 


where  r((j>)  is  a monotonic  increasing  function  with  r(0)  = 0 and  r(l)  = 1;  suitable 
choices  include  r(<f>)  = 0 or  r{<j>)  — 02(3-20).  This  assumption,  in  which  the  density 
does  not  depend  on  temperature  or  pressure,  is  known  as  quasi-incompressibility,  as 
it  still  allows  a nonzero  divergence  of  the  velocity  vector.  This  assumption  places 
a constraint  on  the  form  of  the  underlying  thermodynamic  potentials  [45]  which 
requires  the  underlying  Gibbs  free  energy  (per  unit  mass)  to  have  the  form 

g(T,p,<j))  = go{T,<t>)+{^^,  (19) 

P\4>) 
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where  po  is  a reference  pressure.  Here  we  assume  that  the  function  go(T,  <f>)  has  the 
form 


9o(T,(j))  = 


e0  - cTm  - r(0)L  - - — H{(p) 
4 as 


1-^)-cTln(^)  + ^w’ 


(20) 

in  which  case  the  corresponding  expressions  for  the  internal  energy  and  entropy 
densities  are 


e = e0  + c(T  - TM)  - r(<f>)L  + 

AaE  p[(p) 


s = 


Tm 


e0  - r((f>)L  - 

4 as 


cln 


Tm 


(21) 

(22) 


where  1 /aE  = 1/a— 1/as-  Here  1/a  is  the  height  of  the  double  well  of  the  Gibbs  free 
energy  density  at  T = Tm,  and  l/aE  and  1 /as  are  the  heights  of  the  double  wells 
in  the  internal  energy  and  entropy  densities,  respectively.  The  double  well  potential 
H{(p)  is  a prescribed  function  of  <j>  (see  [56]).  The  quantity  eo  is  a constant  reference 
energy  and  both  the  heat  capacity  per  unit  mass  c and  the  latent  heat  per  unit  mass 
L are  assumed  to  be  constant.  Tm  is  the  melting  point  at  the  reference  pressure 
Po- 

The  tensors  m and  ms  are  given  by 


m = 


ms  = 


-P  + 


4cnr2 


r2(V4>) 


Te2 


_p+_sr2(V0) 


I - e^(r)r(V0)£  <g>  + r, 

I — Te|r(V0)|* 0 V0  + r, 


where  r is  the  viscous  stress  tensor, 
t = n 


T 2 


Vu+  (Vu)J  - -(V  -u)  I 

O 


and  \i  is  the  viscosity,  which  is  a function  of  0, 

p(<t>)  = Psr((p)  + Pl[  1 - r{(t>)}. 


(23) 

(24) 


(25) 


(26) 


By  examining  the  isothermal  one-dimensional  solution  of  the  governing  equa- 
tions at  the  melting  temperature  Tm  with  ps  = Pl,  it  may  be  shown  that  the 
surface  tension  7,  interface  thickness  /,  and  interface  attachment  coefficient  po  are 
related  to  the  phase-field  parameters  by 


7(n)  = 


^dZkirta) 


l(ii)  = eF(TM)T(n) 


Po (n)  = 


6pz,L/r(n) 

TmM 


(27) 

We  will  henceforth  confine  our  attention  to  case  of  isotropic  surface  energies  and 
set  r(V0)  = |V0|;  in  this  case  we  note  that  T(n)  = 1.  It  is  also  convenient  to  define 
an  associated  capillary  length  by  lc  = TmI /{pl[L2 /c]). 
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2.2  Dimensionless  Governing  Equations 


We  non-dimensionalize  the  governing  equations  by  introducing  the  following  di- 
mensionless variables,  which  we  denote  with  a prime: 


I — IqOC  , 


ll  , 
t = — t'. 
KL 


u — 


K-L 

loU' 


m = 


PL^j 

ll 


m’,  p = p0+  —pr^P 
Iq 


(28) 


T = TM  + —O',  p = plp',  P — Plp'  , k = kLk'.  (29) 

c 

Here  the  reference  length  scale  l0  is  a typical  length  scale  associated  with  the  in- 
terface shape,  such  as  a dendrite  tip  radius,  the  reference  time  scale  is  Iq/kl,  and 
the  reference  velocity  scale  is  U = kl/Iq,  where  kl  is  the  thermal  diffusivity  in  the 
bulk  liquid  phase.  The  dimensionless  governing  equations  are 


!+V.(/>u)=0, 


(30) 


Du 

p-=—  = V • m, 
y Dt 


(31) 


e2M^y  = e2(l  + ad)X72(j)  — p 


-(l+WH'W+Xer'W  + ^p—  - 

2 7 <90  \pj  J 


,(32) 


pTt=v'  + e2^  + U'  (33) 

where,  for  simplicity,  we  have  omitted  the  primes  on  the  dimensionless  variables. 
The  dimensionless  internal  energy  density  is  given  by 

e — 0 — r{(t>)  + ^P(0)  — -p— (34) 

2 A7  p 

and  the  dimensionless  stress  tensors  are 

m = crv  4-  (1  + aO)^  + Prr,  (35) 

ms  = crp  + a(9  + S-1)^  + Prr,  (36) 


where  Pr  — ul/ kz.  is  the  Prandtl  number  of  the  liquid  phase,  is  the  kinematic 
viscosity  of  the  liquid  phase  and 


crv  - - {p  + P*)  I , 

1 


<7^  — 76 


r = p(0) 


|V0|2I- V0 


Vu  + (Vu 


-7 


The  source  term  in  the  energy  equation  is 

H = : Vu  = [-  (p  + P*)  I + a 

A7  A7  L 


(37) 

<8>  V0 

(38) 

- |(V -ff)  / . 

(39) 

(0  -1-  S-1)  0^  -1-  Prr]  : Vu, 

(40) 
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and  the  dimensionless  parameters  are  given  by 

aL  „ 1 


a = 


ce2F(TM)' 


P = 


<*scTm  ’ 


6 = 


2 cleL  ’ 


y = 


/2LPl  2aLe2F(T m)  ’ 


(41) 


-1  x__L 

C /o’  6/c’ 


P = 


Po*o 

PL«| 


5 = 


cTm 


(42) 


7 = 


dlplcL2  _ 6/07 
k2lcTm 


Plk2l 


M = 


klMTmc  _ [kl/Ic] 

6 pLllcL2  ~ [ L/c\fi0 ’ 


(43) 


and 


W) 


^(0), 


/i(0)  = l+(  — ~ 1 
\PL 


r(0), 


Ps 


P(0)  = 1+(^“  ~ 1J  r(^)' 
(44) 

We  note  that  the  parameter  7 is  related  to  the  capillary  (or  crispation)  number  Ca 
by  7 = 6 Pr/Ca,  where  Ca  = o7)- 

In  the  absence  of  flow  these  equations  reduce  to  the  generalized  phase-field 
equations  recently  studied  in  Refs.  [52]  and  [53].  The  leading-order  free  boundary 
problem  that  emerges  from  a sharp-interface  limit  of  these  equations  depends  on 
the  distinguished  limit  that  is  taken  [54,55,52].  The  limit  in  which  A = 0(1)  as 
e -)•  0 corresponds  to  the  so-called  ‘thin  interface’  limit  studied  by  Karma  and 
Rappell  [55].  In  this  analysis  when  the  thermal  conductivities  of  the  solid  and 
liquid  phases  are  unequal  the  leading-order  temperature  is  discontinuous  across  the 
interface  and  the  leading-order  modified  Gibbs-Thomson  equation  contains  terms 
dependent  on  the  interfacial  temperature  gradients.  However,  if  A = 0(e)  as  e — > 0, 
the  so-called  ‘classical’  limit,  the  temperature  is  continuous  across  the  interface 
at  leading  order  and  a nonlinear  form  of  the  modified  Gibbs-Thomson  equation  is 
obtained  at  leading  order.  However,  if  we  formally  set  the  coefficients  a,  (3,  5,  and 
v to  zero,  thereby  omitting  the  nonstandard  terms  in  the  generalized  phase-field 
equations,  then  the  classical  sharp  interface  analysis  recovers,  at  leading  order,  a 
standard  free  boundary  problem  in  which  the  interfacial  temperature  is  continuous 
and  the  conventional  modified  Gibbs-Thomson  equation  is  obtained. 

Here  we  study  a simplified  form  of  our  model  by  setting  the  constants  a,  /?, 

6,  and  v to  zero,  and  we  neglect  the  source  term  T-L  in  the  energy  equation.  The 
dimensionless  governing  equations  of  the  simplified  model  are 


dp 

a?  + v 


(pu)  = 0, 


(45) 


Du 

P^  = V-m, 


(46) 


e2M^-  --  e2V20  — p 
Dt 


l-H1  (cfr)  + \0r' ((j>)  + 8 (l 


7 <90  \p)  J 


(47) 


p^  = v-[<3(0,v«]. 


(48) 
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where  the  dimensionless  stress  tensor  is 


m = crp  + + Prr. 


(49) 


We  have  recently  examined  the  full  system  of  governing  equations,  including 
flow,  in  the  sharp-interface  limit  [56].  Our  investigations  reveal  that,  in  the  classical 
sharp-interface  limit,  the  boundary  conditions  at  a sharp  interface  in  equilibrium 
comprise  the  normal  stress  balance  including  surface  tension  and  the  Clausius- 
Clayperon  equation  all  under  isothermal  conditions.  In  the  nonequilibrium  case  we 
find  hydrodynamic  conditions  on  the  normal  and  tangential  velocities  represent- 
ing the  conservation  of  mass  and  the  no-slip  condition.  Jumps  conditions  on  the 
normal  and  tangential  stresses  are  also  obtained.  The  temperature  is  found  to  be 
continuous  across  the  interface  while  the  jump  in  heat  flux  across  the  interface  is 
modified  by  nonequilibrium  effects.  The  temperature  of  the  interface  found  to  obey 
a nonequilibrium  version  of  the  Clausius-Clapeyron  relation. 

3 Model  Computations 

We  now  describe  computations,  based  upon  the  phase-field  model  given  by  equa- 
tions (45),  (46),  (47),  (48)  and  (49),  that  represent  the  density  change  flow  associ- 
ated with  the  growth  of  a dendrite  from  an  undercooled  melt.  To  proceed  we  make 
a number  of  additional  approximations  in  order  to  develop  a simplified  phase-field 
model  that  captures  the  qualitative  features  of  this  situation.  First,  we  consider 
the  dendrite  to  be  two-dimensional  and  growing  in  a uniform  thin  gap  of  width  d 
between  two  thermally  insulated  flat  plates.  This  allows  us  to  ignore  the  effects  of 
inertia  and  to  model  the  flow  using  a Hele-Shaw  approximation.  The  momentum 
equation  may  then  be  written  as 


- eV  • [V0  ® V0]  4-  V • [/x(0)t]  = 0.  (50) 

6 


V 


In  the  absence  of  flow  it  is  known  that  it  is  essential  to  include  surface  energy 
anisotropy  in  order  to  compute  dendritic  structures  using  a phase-field  model  [14]. 
We  will  accordingly  retain  anisotropic  surface  energy  terms  in  the  phase-field  equa- 
tion alone.  Specifically,  an  isotropic  surface  energy  term  is  used  in  the  momentum 
equation  (50)  while  the  phase-field  equation,  given  by  (47),  is  modified  to  allow  for 
anisotropic  surface  energy  by  using  the  Cahn-Hoffman  ^-vector  formalism 


= e2V  • [r(V)|  - p ijT(«  + W(« 


(51) 


so  that  the  direct  effect  of  anisotropy  is  upon  the  interfacial  surface  energy  rather 
than  the  flow.  We  note  we  have  also  omitted  the  pressure  dependence  in  the  free 
energy  term  in  the  phase-field  equation.  This  is  a reasonable  approximation  for 
density-driven  flows,  as  evidenced  by  the  insignificant  variations  of  the  melting 
temperature  due  to  pressure  fluctuations  in  the  Clausius-Clayperon  relation  under 
these  conditions. 

In  order  to  simplify  the  system  further  we  make  the  approximation  V • (V0  <g> 
V0)  « V(|V0|2).  This  approximation  is  exact  in  one  space  dimension,  but  not 
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higher  dimensions.  We  justify  it  by  observing  that  the  phase-field  variable  only 
changes  in  the  thin  interfacial  regions  where  it  depends  primarily  on  the  perpendic- 
ular distance  through  the  layer  and  hence  is  approximately  one  dimensional.  Using 
this  approximation  the  momentum  equation  (50)  becomes 

-Vp+^V-[/z(# r]  = 0,  (52) 

where 

P=h>+5|V<A|2.  (53) 

7 ^ 


We  now  integrate  these  equations  across  the  narrow  gap  of  width  d <£  1 [58]  to 
obtain 


u = 


-3d2  / z\  ( d - z\ 

ca  n(<t>)  \d)  v^d- ; 


(54) 


where  here  and  below  the  operator  V acts  in  the  plane  parallel  to  the  thin  gap.  We 
now  apply  the  continuity  equation  (45)  to  find  that  p satisfies 


V2" 


P = 


d2 


p'W 


p(<t>)  dt  [p(<f>)  p(0). 


V</>  • Vp. 


(55) 


In  our  numerical  computations  we  solved  the  energy  equation  (48),  phase-field 
equation  (51)  and  pressure  equation  (55)  using  VLUGR  [57].  This  freely  distributed 
package  is  designed  to  solve  systems  of  parabolic  partial  differential  equations  in 
which  the  solution  exhibits  regions  in  space  with  large  gradients.  It  employs  a 
finite-difference  discretization  allied  to  local  grid  refinement  and  a variable  time 
step  integration  of  the  underlying  discretized  equations. 

Computations  were  conducted  on  the  rectangular  domain  [0,  X]  x [0,  T].  Neu- 
mann boundary  conditions  were  employed  on  <j>  and  6 on  all  four  sides.  However, 
for  the  pressure,  Neumann  boundary  conditions  were  only  invoked  on  the  sides 
x = 0 and  y = 0,  with  Dirichlet  boundary  conditions  on  the  other  two  sides.  The 
initial  condition  represented  a small  circular  solid  region  centered  on  the  origin  in 
a uniformly  undercooled  melt  with  dimensionless  temperature  T = c[Tm  — To)/ L, 
where  Tq  is  the  initial  dimensional  temperature. 

The  governing  equations  were  solved  with  p{4>)  given  by  Eqn.  (44),  r'(0)  = 
3O02(1  — (p)2  and  H(cj))  = <j)2(  1 — d>)2.  The  surface  energy  had  a four-fold  anisotropy 
with  T(n)  = 1 + 0.005  cos(40),  where  n = V0/|V0|  is  a unit  vector  in  the  ( x,y ) 
plane  and  © is  the  angle  between  n and  the  x-axis.  The  values  of  the  dimensionless 
parameters  used  in  the  computations  are  given  by  Ca  = 30,  p-s/pl  — 1,  ks/ki  = 1, 
A = 7.5,  and  M = 10. 

In  Figure  1 we  display  the  results  of  a typical  computation  in  which  ps/ Pl  = 0.9, 
X = 1,  and  Y = 3.  This  figure  shows  the  pressure  field,  the  velocity  field  and  the 
phase  field  at  time  t = 0.3.  The  solid  curves  are  isobars,  the  arrows  represent  the 
local  velocity,  and  the  shading  indicates  the  phase  field.  The  x and  y axes  represent 
planes  of  symmetry  in  the  calculation,  although  since  1^7  the  resulting  shape 
has  two- fold  but  not  four- fold  symmetry  due  to  the  presence  of  the  sidewalls.  The 
dendrite  growing  in  the  x direction  has  a blunter  tip  than  the  dendrite  growing  in 


Interfaces  for  the  Twenty-First  Century 


10 


Figure  1.  The  phase-field,  pressure  field  and  velocity  field  for  a computation  at  £ = 0.3  with 
Ps/Pl  — 0.9  on  a 1 x 3 domain  with  four-fold  anisotropy.  The  colors  indicate  the  liquid  (</>  = 0) 
and  solid  (<p  — 1}  region,  the  solid  curves  are  the  isobars,  and  the  small  arrows  represent  the 
velocity  field. 


the  y direction  due  to  its  closer  proximity  to  the  sidewall,  which  has  a significant 
effect  on  the  growth  dynamics  at  this  stage. 

Since  the  density  of  the  solid  is  less  than  that  of  the  liquid,  a given  amount  of 
material  will  expand  upon  solidification,  which  drives  a flow  away  from  the  interface 
into  the  melt.  For  a sharp-interface  model,  the  conservation  of  mass  boundary 
condition  takes  the  form 

un  = ~vn  ~ 1^  , (56) 

where  un  = n • u is  the  normal  component  of  the  fluid  velocity  at  the  interface,  and 
vn  is  the  normal  velocity  of  the  interface.  Thus  for  solidification  with  vn  > 0,  the 
flow  is  away  from  the  interface  ( un  > 0)  for  ps  < Pl,  and  is  toward  the  interface 
(un  < 0)  for  ps  > Pl-  The  computation  shows  that  the  flow  is  greatest  in  the 
vicinity  of  the  tip  of  the  dendrite.  For  this  calculation  with  equal  viscosities  in  the 
solid  and  liquid  phases,  there  is  also  a significant  flow  in  the  solid  region.  This 
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artifact  is  reduced  for  computations  with  ps/pl  1;  here  we  are  illustrating  an 
extreme  example  of  this  effect. 


Figure  2.  The  phase-field,  pressure  field  and  velocity  field  at  t — 0.75  for  a computation  with 
Ps/pL  = 1.1  on  a 2 x 4 domain  with  four-fold  anisotropy.  The  colors  indicate  the  liquid  (<£  = 0) 
and  solid  {<f>  — 1)  regions,  the  solid  curves  are  the  isobars,  and  the  small  arrows  represent  the 
velocity  field. 


Figure  2 shows  a similar  situation  but  with  the  density  in  the  solid  greater 
than  that  of  the  liquid,  Ps/pl  — 1.1,  with  X = 2 and  Y = 4 at  a time  t = 0.75. 
In  this  case  the  advection  is  toward  the  interface,  as  expected.  Between  the  two 
dendrite  tips  is  a narrow  liquid  intrusion  where  little  solidification  is  taking  place; 
the  flow  velocities  that  are  induced  by  the  density  change  upon  solidification  are 
correspondingly  small  in  this  region. 

4 Conclusions 

In  this  paper  we  have  shown  that  computations  based  on  a simplified  form  of  a 
recent  phase-field  model  that  includes  convection  [1]  exhibits  numerical  solutions 
which  show  the  expected  physical  behavior.  In  particular,  we  considered  the  growth 
of  a dendrite  in  a thin  gap  between  two  thermally  insulated  plates  and  allowed  the 
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density  of  the  solid  and  liquid  phases  to  be  different.  We  found  that  the  flow  was 
directed  towards  or  away  from  the  dendrite  depending  on  whether  the  density  of 
the  solid  phase  was  greater  or  less  than  that  of  the  liquid  phase,  respectively. 

Our  model  for  solidification  with  convection  is  derived  using  the  formalism  of  ir- 
reversible thermodynamics,  and  allows  the  systematic  incorporation  of  a consistent 
thermodynamic  description  of  the  two-phase  system.  It  allows  a unified  treat- 
ment of  both  equilibrium  and  non-equilibrium  effects  in  a single  set  of  governing 
equations.  Sharp-interface  limits  of  the  diffuse-interface  description  then  lead  to 
boundary  conditions  of  the  solid-liquid  interface  which  recover  the  usual  conditions 
at  equilibrium,  and  provide  thermodynamically-consistent  generalizations  of  these 
conditions  under  non-equilibrium  conditions  [56].  The  detailed  nature  of  the  non- 
equilibrium conditions  at  the  interface  can  be  sensitive  to  the  specific  forms  that 
are  assumed  to  describe  the  variation  of  the  thermophysical  parameters  through 
the  interfacial  region;  for  example,  the  non-equilibrium  solute  trapping  behavior  of 
a diffuse- interface  model  of  a binary  alloy  depends  quantitatively  on  the  exact  form 
that  is  assumed  for  the  variation  of  solute  diffusivity,  near  the  interface  [59]. 

In  this  model  the  solid  is  treated  as  a liquid  with  high  viscosity.  This  allows 
residual  convection  to  occur  in  the  solid,  with  a magnitude  that  is  determined  by 
the  viscosity  ratio.  The  consideration  of  extreme  viscosity  ratios  tends  to  elimi- 
nate velocity  gradients  in  the  solid,  but  allows  states  of  uniform  convection  that 
correspond  to  rigid  body  motion.  This  is  an  attractive  feature  for  dealing  with 
such  issues  as  fragmentation  and  subsequent  motion  of  sidearms  through  the  melt. 
Such  transported  fragments  can  serve  as  sites  for  the  growth  of  independent  grains 
when  the  fragments  are  incorporated  into  the  growing  phase,  which  is  a problem  of 
considerable  technological  importance.  It  is  thus  beneficial  to  have  a model  which 
allows  for  both  topological  changes  in  the  interface  as  well  as  possible  rigid-body 
motion  of  the  solid  phase. 
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